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Abstract 

The evolution of a mixed quantum-classical system is expressed in the mapping formalism where 
discrete quantum states are mapped onto oscillator states, resulting in a phase space description 
of the quantum degrees of freedom. By defining projection operators onto the mapping states 
corresponding to the physical quantum states, it is shown that the mapping quantum-classical 
Liouville operator commutes with the projection operator so that the dynamics is confined to the 
physical space. It is also shown that a trajectory-based solution of this equation can be constructed 
that requires the simulation of an ensemble of entangled trajectories. An approximation to this 
evolution equation which retains only the Poisson bracket contribution to the evolution operator 
does admit a solution in an ensemble of independent trajectories but it is shown that this operator 
does not commute with the projection operators and the dynamics may take the system outside 
the physical space. The dynamical instabilities, utility and domain of validity of this approximate 
dynamics are discussed. The effects are illustrated by simulations on several quantum systems. 

PACS numbers: 
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I. INTRODUCTION 



Since phenomena such as electron and proton transfer dynamics^, excited state relax- 
ation processes^ and energy transport in light harvesting systems^ are quantum in nature, 
the development of theoretical descriptions and simulation methods for quantum many-body 
systems is a central topic of research. Although various techniques can be used to study 
such problems, quantum-classical methods^— , where certain degrees of freedom are singled 
out for a full quantum treatment while other environmental variables are treated classically, 
permit one to investigate large and complex systems that cannot be studied by other means. 

In this article we consider descriptions of the dynamics based on the quantum-classical 
Liouville equation^ (QCLE) and, in particular, its representation in the mapping basis^r— . 
The mapping formalism provides an exact mapping of discrete quantum states onto con- 
tinuous variables^ and in quantum-classical systems leads to phase-space-like evolution 
equations for both quantum and classical degrees of freedom. The mapping basis has been 
used in a number of different quantum-classical formulations, often based on semi-classical 
path integral expressions for the dynamics^ - — . The representation of the quantum-classical 
Liouville equation in the mapping basis leads to an equation of motion whose Liouvillian 
consists of a Poisson bracket term in the full quantum subsystem-classical bath phase space, 
and a more complex term involving second derivatives of quantum phase space variables and 
first derivatives with respect to bath momenta 12 . This latter contribution has been shown 
to be an excess coupling term related to a portion of the back reaction of the quantum 
subsystem on the batbM. 

Various aspects of the QCLE in the mapping basis and properties of its full and approxi- 
mate solutions are discussed in this paper. The solutions of the quantum-classical Liouville 
equation cannot be obtained from the dynamics of an ensemble of independent classical-like 
trajectories^. In the adiabatic basis this equation admits a solution in terms of surface- 
hopping trajectories^— , but other schemes have been used to simulate the dynamics^— . 
When it is expressed in the mapping basis, we show that a solution can be obtained in 
terms of an ensemble of entangled trajectories. The excess coupling gives rise to correlations 
between the dynamics of the quantum mapping degrees of freedom and the bath phase space 
variables that are responsible for the entanglement of the trajectories in the ensemble. The 
derivation of the entangled trajectory picture is similar to that for trajectory solutions of 
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the Wigner-Liouville equation^^. 

If the excess coupling term is dropped and only the Poisson bracket part of the Liouvillian 
is retained, a very simple equation of motion that admits a solution in terms of characteris- 
tics is obtained. Consequently, its solutions can be obtained from simulations of an ensemble 
of independent trajectories evolving under Newtonian dynamics. The set of ordinary dif- 
ferential equations has appeared earlier in mapping formulations based on semi-classical 
path integral formulations of the dynamics^ 1 ^ 1 ^, indicating a close connection between this 
approximation to the quantum-classical Liouville equation and those formulations. The 
solutions of this Poisson bracket approximation to the QCLE, as well as those of other semi- 
classical schemes that use this set of evolution equations, often provide a quantitatively 
accurate description of the dynamics^ 1 ^ 1 ^. However for some systems the solutions are not 
without artifacts and difficulties. Some of these difficulties can be traced to the fact that the 
independent-ensemble dynamics can take the system out of the physical space and inverted 
potentials can appear in the evolution equations, which may lead to instabilities^^^. 

The main results of this paper are as follows: We present derivations of expressions for 
mapping quantum-classical (MQCL) evolution equations and expectation values of operators 
that explicitly show how projection operators onto the physical mapping eigenstates enter 
the formulation. We demonstrate that the MQCL operator commutes with this projection 
operator so that dynamics under this evolution is confined to the physical space. This 
full quantum-classical dynamics in the mapping basis can be simulated by an ensemble of 
entangled trajectories. We also show that when the excess coupling term is neglected the 
resulting Poisson bracket operator no longer commutes with the projection operator so that 
this approximate dynamics can take the system out of the physical space. Given this context, 
we revisit the issue of instabilities in the dynamics of the Poisson bracket approximation 
and discuss the conditions under which such instabilities are likely to arise and lead to 
inaccuracies in the solutions. 

In Sec. [Til we outline the representation of the quantum-classical Liouville equation in the 
mapping basis and show how average values of time dependent observables may be com- 
puted. We also define a projection operator onto the mapping states and show how this 
projector enters the expressions for the expectation values and evolution equations. Sec- 
tion IIHI briefly describes the entangled trajectory solution to the QCLE in the mapping 
basis. This section also shows that when the excess coupling term is neglected, a solution in 
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terms of an ensemble of independent trajectories is possible. In Sec. lIVI the approximate evo- 
lution equation obtained by retaining only the Poisson bracket term in the Liouville operator 
is considered and the dynamical instabilities that can arise in the course of the evolution are 
highlighted. Various aspects of the theoretical analysis that concern the approximate solu- 
tions and resulting instabilities are illustrated by simulations of a number of model systems. 
A brief summary of the main results of the study, along with comments, are given in Sec. |V] 
The Appendices provide material to support the text. In particular, we describe an efficient 
simulation algorithm for the ordinary differential equations that underlie the solutions of 
the Poisson bracket approximation to the QCLE. 

II. QUANTUM-CLASSICAL LIOUVILLE EQUATION: MAPPING, PROJEC- 
TORS AND EXPECTATION VALUES 

The quantum-classical Liouville equation (QCLE), 

^p w (X,t) = -iCp w (X,t), (1) 

describes the time evolution of the density matrix pw(X, t), which is a quantum operator that 
depends on the classical phase space variables X = (R, P) = R 2 , Rn £ , Pi, Pi, Pn e ) 
of the environment. The quantum-classical Liouville operator is defined by 

it = l -[H w , •] - \{{H W , •} - {•, H w }), (2) 

where Hw{X) is the partial Wigner transform of the total Hamiltonian of the system, [•, •] is 
the commutator and {•, •} is the Poisson bracket in the phase space of the classical variables 
X. The total Hamiltonian may be written as the sum of environmental (bath), subsystem 
and coupling terms, H W {X) = H e (X) + h a + V C (R), where H e (X) = P 2 /2M + V e {R) is the 
bath Hamiltonian with V e (R) the bath potential energy, h s = p 2 /2m + V s is the subsystem 
Hamiltonian with p and V s the subsystem momentum and potential energy operators, and 
V C (R) is the coupling potential energy operator. Here m and M are the masses of the 
subsystem and bath particles, respectively. 

The QCLE may be written in the basis, {|A); A = 1, . . . , N}, that spans the quantum 
subsystem space with eigenfunctions defined by the eigenvalue problem, h s \X) = €\\\). 
Taking matrix elements of Eq. we obtain 

^'(X,t) = -iC xx ,, vv ,p${X,t). (3) 
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The Einstein summation convention is used here and in subsequent equations although, on 
occasion, sums will be explicitly written for purposes of clarity. The QCL operator in the 
subsystem basis is^ 

( P 9 n/n. 9 \ S S 

+ {M-dR +F ^ R) -dp) 5 ^' 
1 (. dV c Xu dV»' x '\ d 

where u\\> = (e\ — e\>)/h and F e (R) = —dV e /dR is the force due to molecules in the 
environment. 

The evolution equation for an observable Bw{X), analogous to Eq. (CQ), is 

j t B w (X,t)=i£B w (X,t), (5) 

and its representation in the subsystem basis is analogous to Eq. fl3]) with a change in sign 
on the right side. 



A. Representation in Mapping Basis and Projection Operators 

In the mapping basis^ 1 ^ the |A) eigenf unctions of an iV-state quantum subsystem can be 
replaced with eigenfunctions of iV fictitious harmonic oscillators, \m\), having occupation 
numbers which are limited to or 1: |A) — > \m\) = |0i, ■ • • , 1a, • ■ ■ On)- Creation and anni- 
hilation operators on these states, a\ and a\, respectively, may be defined. For any operator 
B\y{X) whose matrix elements in the subsystem basis are Byy'(X), we may associate a 
mapping basis operator B W (X) — > B m (X), where 

B m (X)=B$'(X)a[a x , (6) 

It is then evident that the matrix element B^' (X) = (\\Bw(X)\\') = (m\\B m (X)\m\') . 

The expression for B^' (X) may also be written in terms of the Wigner transforms in the 
space of the mapping variables. Inserting complete sets of coordinate states {\q), \q')}, and 
making the usual coordinate transformations appropriate for Wigner transforms, (q, q') — > 
(r — z/2, r + z/2), we obtain 

B$\X) = (m x \B m (X)\m y ) = (7) 

/z z z z 

drdz (m\\r - -)(r - -\B m (X)\r + -)(r+ -\my) 
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Another form for the matrix element can be obtained by inserting the Wigner transform of 
an operator and its inverse as 

(r - Z -\B m {X)\r + |> = — Ip J dp e-^ h B m (X), 

B m (X) = J dz e*-*/ R (r - ||5 m (X)|r + ~). (8) 

Here X = (x, X) are the extended phase space coordinates for the subsystem mapping 
variables, x = (r,p) = (ri, r^,Pi, ■ ■■,Pn), an d the environment, X = (R,P). Making 
these substitutions in Eq. ([7]) we obtain, 

B$\X) = j dx B m (X)g xx ,(x), (9) 



where we have defined- 



1 / z z 

9xx '^ = (2ixK) N J dZ e ~ iP ' Z/h ( r + ^K'X^lr ~ 2^ 

' dze ip - z/h (r-^\m y )(m x \r + J). (10) 



(2Trh) N J v T /v 1 2 

Evaluating the integral we obtain an explicit expression for g xx >(x): 



gM*) = <P(?) (11) 
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r\r x > + PaPa' - i \ r xP\' - ^a'Pa) - -r X \> 



where (j>(x) = {nh)~ N exp (—x 2 /h) is a normalized Gaussian function. Here x 2 = r x r x + p x p x 
in the Einstein summation convention. 

The expression for B m (X) in Eq. flE]) can be simplified by evaluating the integral in the 
Wigner transform. Using the definition of B m (X) in Eq. (jBJ), Eq. (jSJ) may be written as 

B m (X) = B^'(X) J dz e^/\r - ^\a{a x ,\r + |). (12) 

Noting that the factor multiplying B^' (X) is the Wigner transform of a\a X i, (d\a X t)w{x) = 
c X y(x), whose explicit value is 

c XX '{x) = ^[ r A^A' +P\P\> + i(r x p x > - r y p x ) - hS xy ], 

(13) 

we find 

B m {X) = B$\X)c xx ,{x). (14) 
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We may deduce a number of other relations given the definitions stated above. A mapping 
operator B m (X) acts on mapping functions \m\). In this space we have the completeness 
relations V = Yl\=i \ m x)( m \\ — 1> where V is projector onto the complete set of mapping 
states. 40 - Thus, a mapping operator can be written using this projector as 

Bl{X)=VB m {X)V = \m x }(m x \B m (X)\my}(m y \ 

= \m x )B%'(X)(m y \, (15) 

where in the second line we used the equivalence between matrix elements in the subsystem 
and mapping representations given in Eq. (jTJ). We can make use of the Wigner transforms 
defined in Eq. (jSJ) to write these relations in other forms. Using the first equality in Eq. ( TT5|) 

we have 

B*(X) = Jdz e W»( r _ Z -\Bl(X)\r + |> = 

dz e ivz/K( T _ ^\ mx )( mx \ Bm (X)\m x ,}(m x ,\r + |> 

= (2irh) N g x , x (x) J dx' g xy (x')B m (x' , X), 
= VB m (X), (16) 

where we used Eqs. (GO) and The last line defines the projection operator V that projects 
any function of the mapping phase space coordinates, f(x), onto the mapping states, 



Vf(x) = (27cK) N g x , x (x) J dx' g X y(x')f(x'). (17) 

One may verify that V 2 = V since 

(2irh) N J dx g xx >(x)g v > u (x) = 6 Xv 5 x > u >. (18) 

An equivalent expression for B^(X) can be obtained by starting with the last equality 
in Eq. (fT5|) and taking Wigner transforms to find 

B%(X) = jdz e**t*(r - ~\m x )B$' (X)(m x ,\r + *-) 

= (2nh) N g x , x (x)B^(X). (19) 

This result also follows from Eq. (Tl6l) by substituting Eq. ( I14p for B m (X) and using the fact 
that 

dx g xx <{x)c uu ,(x) = 5 X J X > U >. (20) 



Finally, in view of the definition of the projection operator V, in place of Eq. (Q we may 
write 

B%'{X) = J dx Bl{X)g xx ,{x). (21) 

An analogous set of relations apply to the matrix elements of the density operator, 
p$(X) = (X\p w (X)\X') = (m x \pm(X)\m X '), where p m (X) = p^'(X)a[d y . Taking the 
Wigner transform of p m (X) we find 

Pm{Z) = (2^F / dz ^ Z/h ^-l\Prn{X)\T+ Z -) 

' p A /(X)c AV (x). (22) 



(2^)^ 

Likewise, starting from the expression for the projected density. 



pZ(X) = \m x ){mx\pm{X)\m X '){mx'\ 

= \m x )p#\X){m x ,\, (23) 



its Wigner transform is 



f&W = J dpe^l\r- Z -\pl{X)\r+ Z -) 

= V Pm {X), (24) 

which, repeating the steps that gave Eq. (Tl9i) . yields 

ffl = #(%*(i). (25) 
Following the analysis given above that led to Eq. (Q for an operator, and using the relation 

(r - Z -\p m {X)\r + Z -) = j dp e^ h p m (X), (26) 

the evaluation of p^' (X) = (m A |p m (X)|m A /) leads to 

p^\X) = (27rh) N J dx gxx>(x)p m (X) 



(2nh) N J dx gxxix)pl(X). (27) 



These relations allow one to transform operators expressed in the subsystem basis to 
Wigner representations of operators in the basis of mapping states. The projected forms 
of the mapping operators and densities confine these quantities to the physical space and 
this feature plays an important role in the discussions of the nature of dynamics using the 
mapping basis. We now show how these relations enter the expressions for expectation 
values and evolution equations. 
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B. Forms of Operators in the Mapping Subspace 



We first consider the equivalent forms that operators take, provided they are confined to 
the physical mapping space. Since 

(m x \ s ^ala v \m X ') = (m x \m x >), (28) 

Y2u wja v is an identity operator in the mapping space. (Here we include the explicit sum- 
mation on mapping states for clarity.) Using the definition of g XX '{x) in Eq. fllOp . we may 
write the right side of Eq. f[2"8"j) as 

(m x \m y ) = J dx gxy(x). (29) 

The left side of may be evaluated by inserting complete sets of coordinate states and taking 
Wigner transforms so that an equivalent form for Eq. (|28|) is 

/«fc 9 „(,)E<w(*)=/^ (30) 

Thus, we see that 

E c -( x ) = ^ X>2 + pt - h ) = !> ( 31 ) 

V V 

provided it lies inside the gw(x) integral. 

This result has implications for the form of operators in the mapping basis. The matrix 
elements of an operator Bw(X) in the subsystem basis may always be written as a sum of 
trace and traceless contributions, 

B^'(X) = 8 xx ,(TrB w )/N + Bw(X), (32) 

\\> 

where B w (X) is traceless. Inserting this expression into Eq. fl!4p for B m (X), we obtain 

\\i 

B m (X) = (Tr B w )/N + B w (X)c xx ,(x), (33) 

provided B m (X) appears inside the g X y(x) integral. Note that all subsystem matrix elements 
are of this form in view of Eq. ([9}. Here c XX '(x) = ^ [t x t x > + p x p x > + i(r x p X i — r x >p x )] is the 
traceless form of c xx r(x). 

As a special case of these results, we can write the mapping Hamiltonian, H m (X) = 
H$ (X)c XX i(x) in a convenient form. The Hamiltonian matrix elements are given by 

H$'(X) = H e (X)6 xx ,+e x 6 xx , + V c xx '(R) 

= H e (X)5 xx , + h xx '(R), (34) 



which can be written as a sum of trace and traceless contributions, 

H#\X) = (H e {X) + {i:vh)/Ny xxl + h XX \R) 

= H (X)8 XX , + h (R). (35) 

The Hamiltonian H can be written as H = P 2 /2M + Vq(R). From this form for H^' , it 
follows that 

P 2 1 -A A' 

H m (X) = 2M+Vo(R) + ^h (R)(r x r y +p x py), (36) 

\\i 

again, when it appears inside integrals with g\x(x). We have used the fact that h is 
symmetric to simplify the expression for c\\> (x) in this expression. This form of the mapping 
Hamiltonian will play a role in the subsequent discussion. 



C. Expectation values 

Our interest is in the computation of average values of observables, such as electronic 
state populations or coherence, as a function of time. The expression for the expectation 
value of a general observable B W (X) is 



B(t) = dXTr (B w (X)pw(X,t)) 



(37) 



dX B^\X)p x w \X,t) = / dX B#'(X,t)f%f(X) 



,A'A 



j XX' 



A' A/ 



where the trace is taken in the quantum subsystem space. In the last line the time de- 
pendence has been moved from the density matrix to the operator, which also satisfies the 
QCLE. 

The expression for the expectation value can be written in the mapping basis using the 
results in the previous subsection. For example, using Eq. Qj and the first line of Eq. (T5j 
we find 



B(t) 



dX 



dx B m (X,t)gxx'(x) 
(2nH) N J dx' g yx (x')p m (x',X) 
j dX B m (X,t)pl(X) = j dX Bl(X,t)p m (X), 



(38) 



where we have made use of the definition of the projection operator in Eq. (fTTI) in writing 
the second equality. The projection operator can instead be applied to the observable in 
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view of the symmetry in the expression and the resulting form is given in the last equality. 
We may write other equivalent forms for the expectation value. Starting from the second 
equality in Eq. ( 157)) involving the time evolved density and the time independent operator, 
we obtain 



From a computational point of view, the penultimate equality in Eq. (155)) is most convenient 
since its evaluation entails sampling from the initial value of the projected density and time 
evolution of the operator. 

D. Equations of motion 

The most convenient form of the expectation value requires a knowledge of B m (X,t) = 
B$ (X, t)cw(x). Of course, if the solution to the QCLE in the subsystem basis, B$ (X, t), 
is known, this definition can be used directly to construct B m (X,t); however, the utility 
of the mapping basis representation lies in the fact that one can construct and solve the 
equation of motion for B m (X,t) directly. The derivation of the evolution equation was 
given earlier .— Here, we derive the evolution equations by taking account of the properties of 
mapping operators under integrals of gw(x) in order to make connection with the projected 
forms of operators and densities. This will allow us to explore the domain of validity of the 
resulting equations. 

The QCLE for an observable is expressed in the subsystem basis by taking matrix ele- 
ments of the abstract equation dBw(t)/dt = iCBy/if) with iC defined in Eq. (j2j): 




(39) 



j t (X\B w (X,t)\\') = -±(\\[H w ,B w (X,t)]\X) 
+^(\\({H w ,B w (X,t)} - {B w (X,t),H w })\\'). 



(40) 



We may write this equation in terms of mapping variables using Eq. ()9j) as 




(41) 
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The mapping variables occur inside integrals of g X y(x) integral; i.e., they are projected 
onto the space of mapping states. Since the commutator and Poisson bracket terms in this 
equation involve products of operators, we must obtain the mapping form for a product of 
operators Aw(X)Bw(X). The most direct way to make this transformation is to consider 
the product of operators as they appear in the subsystem basis and then use Eq. §§§ for 
each matrix element: 

A%(X)B$(X) = J dxA m (x,X)g Xu (x) (42) 
x / dx' g vX .{x')B m (x',X). 



This expression does not lead to a useful form for the equations of motion. Instead we may 
write 

A%(X)B#\X) = (X\A W (X)B W (X)\X') 

= (m x \A m (X)B m (X)\m x ,) (43) 
dx g xx/ (x)(A m (X)B m (X)) w (X) 

Given that the Wigner transform of a product of operators is 

(A m (X)B m (X)) w = A m (x,X)e hA -/ 2t B m (x,X), (44) 

where A m = v p ■ V r — ■ V p is the negative of the Poisson bracket operator on the mapping 
phase space coordinates, we obtain 

A%(X)B#\X) = (45) 
dx g xx ,(x) (A m (X)e hAm/2l B m (X) 



In Appendix A we establish the equality between this form for the matrix product and that 
given in Eq. (T42l . Inserting this result into Eq. fj4T|) . expanding the exponential operator 
and noting that the mapping Hamiltonian is a quadratic function of the mapping phase 



space coordinates, we obtain (details of the derivation are given in Ref. [121 ]) 

/d 
dx g xx ,(x)(-B m {X,t) = i£ m B m {X,t)), (46) 

where the mapping quantum-classical Liouville (MQCL) operator is given by the sum of two 
contributions: 

zC m = iC p m B + iC' m . (47) 
12 



The Liouville operator iC^ B has a Poisson bracket form, 

nPB r tt ^ ^ ( 9 d 

<9# m 9 P 9 \ 



<9P <9P M dRJ' ^ 
where {-,-}x denotes a Poisson bracket in the full mapping-environment phase space of the 
system, while 

_h dh xx ' / d 2 d 2 \ _d_ 

In writing this form of the mapping Liouville operator we used the expression for the Hamil- 
tonian given in Eq. (I36I) . This is allowed since by Eq. ( 142|) the operators appear inside gw 
integrals. 

The formal solution of the equation of motion for B m (X,t) is B m (X,t) = e tCmt B m (X). 
The expectation value of this operator is given by (see Eq. (138]) ) 

W) = JdX (e^BU^pKX) = (50) 

j dX B m (X)e- iCmt pZ (X) = J dX B m (X)pl(X,t), 

where the evolution operator has been moved to act on the projected density using integra- 
tion by parts. Thus, we see that the projected density satisfies 

d 

foPm(X,t) = -i£mPm(t)- (51) 

Making use of the above results, we can establish relations among the various forms of 
the expectation values and the dynamics projected onto the physical mapping states. From 
Eqs. (J39D and (EDI) we have the relation f dX B m (X)p v m {X ,t) = j dX B^{X)p m {X, t). 
Differentiating both sides with respect to time and using the MQCLE we may write this 
equality as 

/ dX B m (X)i£ m V Pm (X,t) (52) 
dX B m {X)ViC m p m {X,t). 



This identity, which is confirmed by direct calculation using the explicit form of i£ m in 
Appendix B, shows that iC m commutes with the projection operator. Thus, evolution 
under the MQCL operator is confined to the physical mapping space. 

13 



III. TRAJECTORY DESCRIPTION OF DYNAMICS 



A variety of simulation schemes have been constructed for the solution of the QCLE, some 
involving trajectory based solutions^— . These schemes involve either ensembles of 
surface-hopping trajectories or correlations among the trajectories. A solution in terms of an 
ensemble of independent trajectories evolving by Netwonian-like equations is not possible^. 

A. Ensemble of entangled trajectories 

A trajectory based solution of the MQCLE can also be constructed but the trajectories 
comprising the ensemble are not independent. Such entangled trajectory solutions have been 
discussed by Donoso, Zheng and Martens^ 1 ^ for the Wigner transformed quantum Liouville 
equation. While our starting equation is very different, a similar strategy can be used to 
derive a set of equations of motion for an ensemble of entangled trajectories. 

The MQCLE ([!]) can be written as a continuity equation in the full (mapping plus envi- 
ronment) phase space as 



The second equality in Eq. f )53|) defines the phase space velocity field v(X; p^(X, £)) through 
j(X,t) = v(X; p^(X, t))p^(X, £)), which is a functional of the full phase space density. 
We seek a solution in terms of an ensemble of Af trajectories, p^(X,t) = 

Wj is the initial weight of trajectory i in the ensemble. 
To find the equations of motion for the trajectories, consider the phase space average of the 
product of an arbitrary function f{X) with Eq. (153j) : 




(53) 




(54) 




(55) 



14 



where we have carried out an integration by parts to obtain the right side of the equality. 
Substitution of the ansatz for the phase space density into this equation gives 



E 

i=l 



df(Xi(t)) 
' dXi(t) 



0. 



(56) 



from which it follows that the trajectories satisfy the evolution equations, Xi(t) 
v(Xi(t); Pm{Xi{t)))- More explicitly we have 

dH m . dH m ■ dH m 



P 



dHm 



Px 



+ 



h dh 



XX' 



R 

d 2 



dP 



+ 



(57) 



d 2 



v 

Pm' 



dR 8p^ dR \drydrx ' dpydpx- 
The second term in the environmental momentum equation couples the dynamics of all 
members of the ensemble since it involves the phase space density. 



B. Ensemble of independent trajectories 

If the last term in the P equation is dropped we recover simple Newtonian evolution 
equations: 

(58) 





9H m 


dpx 


dH m 


dt 


dpx ' 


dt 


dr x 


dR 


dH m 


dP 


dH m 


~dT 


dP ' 


dt 


dR ' 



This result also follows from the fact that neglect of the last term in the P equation cor- 
responds to the neglect of the last term in the formula for iC m in Eq. fj47|) . Thus, in this 
approximation 

§£/C{x,t) = {n m , P v m } x = -iC p m B P v m {x,t), (59) 

which we call the Poisson bracket mapping equation (PBME). Since the approximate evolu- 
tion has a Poisson bracket form, it admits a solution in characteristics and the corresponding 
ordinary differential equations are those above in Eq. fl5Bl 12 . 
In contrast to Eq. (|52l) . in Appendix B we show that 

j dX B m {X)iL p m B V Pm {X) (60) 
^ J dX B m {X)ViC p m B Pm {X). 
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Consequently, the Poisson bracket mapping operator iC^ does not commute with the 
projection operator. Therefore, unlike the evolution under the full MQCL operator, the 
evolution prescribed by the PBM operator may take the dynamics out of the physical space. 

As will be seen shortly, one consequence of the dynamics leaving the physically relevant 
regions of phase space is a lack of stability of trajectories due to inversion of the potential 
for bath coordinates. It is therefore important to minimize artificial instabilities arising due 
to the use of too large a time step in numerical methods of solving the evolution equations. 
We note that as in the case of Brownian motion, the bath coordinates typically evolve on 
a much longer time scale than the subsystem phase space coordinates, as can be seen from 
a scaling analysis of the equations of motion Eq. (I58p in terms of the dimensionless mass 
ratio e = (m/M) 1 / 2 . As a consequence, one might expect that the motion of the subsystem 
limits the size of the time step utilized in the integration scheme, and small time steps must 
be chosen to deal with regions of phase space in which rapid changes in population occur. 
In Appendix C we show that an integrator may be designed using the exact solution of the 
subsystem equations of motion when the bath position is held fixed. Using this integrator, 
numerical instabilities are minimized, allowing us to focus on true instabilities inherent in 
the physical system arising from the PBME approximation. 

We also remark that although these equations of motion have been derived from an 
approximation to QCL dynamics in the mapping basis, they also appear in the in the semi- 
classical path integral investigations of quantum dynamics by Stock and Thoss^ 1 ^ and in 
the linearized semiclassical-initial value representation (LSC-IVR) of Miller— These 
results indicate that LSC-IVR dynamics is closely related to this approximate form of the 
QCLE. Connections between QCL dynamics and linearized path integral formulations have 
been discussed in the literature^ 1 ^. The utility of this approximation to the QCLE hinges 
on the form of the Hamiltonian and the manner in which expectation values are computed. 
These issues are also discussed in the next section. 
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IV. DYNAMICAL INSTABILITIES IN APPROXIMATE EVOLUTION EQUA- 
TIONS 



In Sec. Ill Bl we showed that the mapping Hamiltonian, 

H m {X) = H$ '(X)c X y(x) = (61) 
H$'(X) — {r x rx> + P\P\> + i{r\P\> - r x >P\) - h5 X y], 

could be written in the equivalent form given in Eq. ( )36|) . provided the Hamiltonian operator 
appears inside the gw integral; i.e., is projected onto the physical space. In view of Eq. 
and its equivalence to Eq. (145 p . this condition is satisfied for evolution under the MQCLE. 
Evolution under MQCL dynamics is confined to the physical space and the two forms of the 
Hamiltonian will yield equivalent results. In this section we discuss instabilities that may 
arise in approximations to the MQCL as a result of the dynamics taking the system outside 
of the physical space. Problems associated with the lack of confinement to the physical space 
in other mapping formulations have been discussed earlier by Thoss and Stock~. Here we 
reconsider some aspects of these issues in the context of the QCL formulation. 

While different forms of the mapping Hamiltonian are equivalent in the mapping subspace, 
should the dynamics take the system out of this space, the evolution generated by the 
different Hamiltonian forms will not be the same. Indeed, depending on precise form of 
the dynamics, instabilities can arise that depend on the form of the Hamiltonian that is 
employed. In particular, from the structure of H m in Eq. ( 16T1) . one can see that it is possible 
encounter "inverted" potentials if the quantity in square brackets is negative. This problem 
has appeared in approximate schemes based on the mapping formulation and suggestions for 
its partial remedy have been suggested^ 1 ^ 1 ^. Such investigations have led to the observation 
that the form of H m in Eq. (1361) . where the resolution of the identity is used to simplify the 
Hamiltonian form, provides the best results. 

Even if such inverted potentials are not present at the initial phase points of the trajec- 
tories representing the evolution of the density matrix, they may still arise in the course of 
approximate evolution that may take the system outside the physical space; for example, un- 
der PBME dynamics. To investigate the conditions under which unstable dynamics appear, 
consider systems that have localized regions of strong coupling among diabatic states and 
asymptotic regions where such coupling vanishes. The Hamiltonian matrix is approximately 
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diagonal in the asymptotic regions and in such regions H m takes the form, 



p2 



+ v (R) + J2h T 



p2 



(62) 



H, 



A = 



+ v a 



2M 



2M 



x 



where we have defined I\ = -^{r\+p 2 ) ). The second equality defines the effective asymptotic 
potential energy V asy . Since {T\,H m }x = 0, the T x are conserved in the asymptotic regions 
and can be considered constants. 

The effective asymptotic potential energy can be written in the equivalent form, 



where ^ a Ta = T and AT\ = T\ — T/N, which satisfies J2\ = 0. From this equation 
we see that if the matrix element h xx dominates asymptotically, an inverted potential will 
be possible if AI\ < — j?. If instead V dominates asymptotically, no instability will occur. 
Likewise, if another h x ' x ' grows more quickly asymptotically, and does not lead to an inverted 
potential contribution, it will compensate for the inversion due to the h xx term. Note that not 
all h xx terms can give rise to inverted contributions at the same time because Y2\ = 0. 
An interesting case occurs when all h xx grow asymptotically in the same way, e.g. as h. In 
that case, the asymptotic potential takes the form V asy = Vo(-R), which is never inverted. 

Thus, if not all h xx have the same asymptotic behavior and Vq is not asymptotically 
dominant, then it is possible that inverted potentials may occur. In these cases, even if 
the initial condition is such that an inverted potential does not exist, as the system moves 
through the coupling region and into the asymptotic region, one can encounter cases where 
AT\ < — jj, which may result in an inverted effective potential. 

A. Simulations of the dynamics 

While the evolution prescribed by the PBME in Eq. (1591) may take the system outside 
the physical mapping space resulting in dynamical instabilities that could affect the quality 
of the solutions, simulations on a variety of systems has shown that often very accurate 
results can be obtained at a computational cost that is far less than that for simulations 
of the full QCLE. For example, accurate results for the spin-boson system^, simple curve 
crossing models^ and the room temperature excitation transfer in the Fenna-Mathews-Olsen 
light harvesting complex^ have been obtained using this method. In this section we have 




(63) 
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chosen examples to illustrate cases where the simulations of the PBME exhibit more serious 
deviations from the solutions of the full QCLE and exact quantum dynamics as a result of 
the effects discussed above. 

1. Curve crossing dynamics: nuclear momentum distributions 
The simple curve crossing model 4 ^ with Hamiltonian 

H m (X) = P 2 /2M+h XX '(R)(r x r x , +p x px>), 
t 1 = -t 2 = A[l-e- B W]R/\R\ 7 

t 2 =t l = Ce~ DR \ (64) 

is one of the common benchmark cases for quantum dynamics. In this model Hfy is traceless 
so the forms of the mapping Hamiltonian in Eqs. ( 136]) and f ISTj) are identical. Quantitatively 
accurate results for population transfer and coherence have been obtained for this system 
using PBME dynamics 14 , so we focus instead on the properties of the nuclear degrees of 
freedom. 

We have shown 1 - 4 that only part of the back coupling of the quantum subsystem on the 
bath is accounted for in this formulation so that the evolution of the classical degrees of 
freedom may differ from that in the full QCLE. Simulations of this model system 20 using 
LSC-IVR approximations to path integral dynamics have shown that the nuclear momentum 
distribution, after the system passes through the avoided crossing, has single peak. More 
accurate simulations based on the forward-backward (FB)-IVR yield a double-peak structure 
in accord with exact quantum results. As the system passes through the avoided crossing 
and the coupling vanishes, the nuclear momenta have characteristically different values in 
the two asymptotic states giving rise to a bimodal distribution. The single-peaked structure 
of the LSC-IVR simulations was attributed to the mean-field nature of the nuclear dynamics 
in this approximation to the dynamics 20 . 

Here we present comparisons of the nuclear momentum distributions obtained from the 
simulations of the QCLE using a Trotter-based algorithm 30 and its approximation by the 
PBME. We expect the PBME to yield results similar to those of LSC-IVR since the evolution 
equations are similar in these approximations^. The momentum distributions are shown in 

Fig.m 
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The PBME simulations do indeed yield a momentum distribution with a single peak. 
The full QCLE simulations are able to reproduce the correct double-peak structure of this 
distribution^, indicating that the failure of the PBME to capture this effect is due to 
the approximations made to obtain this evolution equation, and not the underlying QCL 
description. 



2. Conical Intersection Model 



A two-level, two-mode quantum model for the coupled vibronic states of a linear ABA 
triatomic molecule has been constructed by Ferretti, Lami, and Villiani (FLV) 52 - 1 ^ in their 
investigation of the dynamics near a conical intersection. The nuclei are described using 
two vibrational degrees of freedom: a symmetric stretch, X, the tuning coordinate and an 
anti-symmetric stretch coupling coordinate, Y. We denote the mapping Hamiltonian for 
this model by H!^(R S , P s , x) whose form is given by Eq. fl36l) with 

H ^ = {£k + £k) + i < 65 > 

+ l -M Y u 2 Y Y 2 + \m x u x \(X - X,f + (X- X 2 f], 



and 



-rll t-22 



1 



2 



1 



h = -h = -M x u x X(X 2 - X x ) + -(Xf - Xi) 



2 



h 12 = h 21 = lY e- a{x ~ x ^e- pY \ (66) 

In these equations R s = (X,Y) while P s = (Px,Py), (M X ,M Y ) and (wj,wy) are the 
momenta, masses and frequencies of the X and Y degrees of freedom.— If the FLV model 
is bilinearly coupled to a bath of independent harmonic oscillators the Hamiltonian has the 
form, 

H m (R,P) = H s m (R s ,P s ,x) (67) 



Nb Pf Mm? 



The coordinates and momenta of each bath oscillator with mass Mj are (Rj, Pj) and and iV^ 
is the number of oscillators. The coupling constants and frequencies, Cj and Uj, correspond 
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with those of a harmonic bath with an Ohmic spectral density. The dynamics of the FLV 



model, with and without coupling to the oscillator bath, was studied in detail Ref. |55[ using 
Trotter-based simulations of the QCLE. Below we present results on this model using the 
approximate PBME dynamics. 

In Fig. [2] we compare the PBME results for ground adiabatic state populations P$o at 
t = 50 fs as a function of the coupling strength, 7, with exact quantum and full QCLE 
results. 

We see that general shape, including the appearance of a minimum and maximum in the 
probability as 7 increases, is captured by all methods. The QCLE results reproduce the 
exact quantum results for the FLV population transfer curve in the low coupling range and 
deviate somewhat for intermediate and high values of the coupling. The PBME results are 
less accurate at low coupling strengths and match the QCLE simulations at high coupling. 
While not quantitatively accurate over the full coupling range, the PBME results capture 
the essential physics in these curves. 

It is interesting to examine statistical features of the ensemble of independent trajectories 
that were used to obtain these results. When the calculation were carried out using the 
original mapping Hamiltonian with form in Eq. f l6Tj) we found that 60% of the ensemble was 
initially on an inverted surface. Furthermore, 61% of trajectories in the ensemble experienced 
an inverted surface at least one time step during the evolution and 50% of the trajectories in 
the ensemble diverged. If instead the Hamiltonian with form in Eq. ( )36l) was used 0.1% of the 
ensemble was initially on an inverted surface, 7% of the ensemble experienced an inverted 
surface at least one time step during the evolution and no trajectories in the ensemble 
diverged. In accord with other investigations, these results indicate the sensitivity of the 
approximate evolution equations to the form of the mapping Hamiltonian. Many of the 
effects arising from instabilities can be ameliorated by first separating the Hamiltonian 
matrix into trace and traceless parts and employing the resolution of the identity. 

Figure [3] compares the Px momentum distributions of the FLV model after passage 
through the conical intersection obtained from simulations of the full QCLE and its PBME 
approximation. The figure also presents results for this momentum distribution when the 
FLV model is coupled to a harmonic bath. The QCLE distribution is much narrower than 
that obtained from the PBME simulations and the peak is shifted to somewhat smaller 
momenta. This trend persists when a larger environment is present but the distributions 
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are in much closer accord. This is consistent with the fact that the PBME does not properly 
account for a portion of the influence of the quantum system on its environment. For a 
larger many-body environment the effect of such back coupling will be smaller. 



3. Collinear Reactive Collision Model 

Finally, we consider a two-level, two-mode quantum model^ for the collinear triatomic 
reaction A + BC — > AB + C. The diabatic states of the system are functions of R = (X,Y), 
where X is the distance between atoms B and C, while Y is the distance between atom A 
and the center of mass of the diatomic BC. The mapping Hamiltonian again has the form 
given by Eq. ( |36|) with 



with h = —h . In these equations (Px,Py) and (Mx,My) are the momenta and iner- 
tial masses corresponding to the BC and A — BC degrees of freedom, respectively. This 
model describes two separate diabatic surfaces, and the off-diagonal diabatic coupling matrix 
elements are constant. 

The QCLE for this model has been simulated in the diabatic basis using the multiple 
spawning molecular dynamics method^ and the results are in quantitative agreement with 
numerically exact quantum dynamics^. This system provides an interesting test case since 
the dynamics can, in principle, explore unphysical regions in the model equations. Diver- 
gences occur where the (diagonal) elements of Hamiltonian are large; i.e., for large negative 
values of X and Y — X/2. While the model allows these negative values, physically, they 
represent distances which should not become negative and the model loses its validity in 




-a(X-X ))2 + c - a (Y-X/2-X ) 
-a(F-X/2-X )\2 , D r -a(X-X ) 



(68) 



and 




(69) 
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these regions. Because the potential is large for large, negative values of these coordinates, 
the nonphysical regions are exponentially suppressed if full quantum or full QCL simulations 
are carried out and physically meaningful results can be obtained with this Hamiltonian. 
This is not the case for the dynamics given by Eqs. (|58p due to the instability from the 
inverted potential, and these approximate evolution equations are much more sensitive to 
the form of the potential. 

In order to ensure that the coordinates of the system do not diverge, the model can 
be altered to avoid nonphysical values of the coordinates. A reasonable adjustment of the 
model that keeps the values of X and Y bounded, even in the approximate PBME, is to 
add a steep confining potential, 

V a (R) = D e ( e -*«(^-*0 + e -MY~x/2-x o) ^ ^ (7Q) 

where z and Xq are parameters. We have chosen the following values: z = 4 and Xq = X$/2. 
By denoting the additional potential as V a , the adjusted Hamiltonian is still of the same 
general form, so none of the formalism needs to be changed. We confirmed that this added 
potential does not substantially change the physical problem 57 . 

In the simulations of the reaction dynamics, the initial wave packet was directed towards 
the reaction region by giving it a non-zero y-momentum. This initial momentum can be 
converted to an excess energy, which is roughly the kinetic energy minus the energy of the 
barrier in the reaction. The results of simulations of the PBME are compared with exact 
quantum and full QCLE simulations in Fig. HI This figure plots the reaction probability 
versus the excess energy. The approximate PBME dynamics fails to capture the peaked 
structure of the reaction probability but does yield probabilities which are qualitatively 
comparable to the exact results. We note, however, that if the model Hamiltonian is not 
supplemented with the confining potential, the approximate mapping dynamics diverges and 
no solution is possible. Neither the exact quantum dynamics nor the full QCLE dynamics 
suffers from this problem. This indicates that if the mapping dynamics is not confined to the 
physical space, the instabilities can probe unphysical regions of models with high probability 
and spoil the results. 
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V. SUMMARY AND COMMENTS 



This investigation of the representation of the quantum-classical Liouville equation in the 
mapping basis led to several results. From considerations of how the equations of motion and 
expectation values involve projectors onto the mapping states corresponding to the physical 
space, it was demonstrated that the QCL operator commutes with the projection operator so 
that the dynamics is confined to the physical space. Further, it was shown that a trajectory- 
based solution of this equation entails the simulation of an ensemble of entangled trajectories. 
The development of suitable algorithms for the simulation of entangled trajectories is a topic 
of current research. 

The PBME approximation to the QCLE is closely related to the equations of motion in 
the LSC-IVR approximation to quantum dynamics^^^ 1 ^. It neglects a portion of the 
back coupling of the quantum subsystem on its environment and does admit a solution in 
terms of an ensemble of independent Newtonian-like trajectories, but the dynamics does not 
commute with the projection operator and, thus, the dynamics may take the system outside 
the physical space. This can lead to unstable trajectories arising from inverted potentials in 
the equations of motion. In addition to initially unstable trajectories, dynamical instabilities 
can arise in the course of the evolution. As in other studies^^ 1 ^, these instabilities are 
partially removed by a judicious choice of mapping Hamiltonian. In this circumstance the 
PBME equation yields qualitatively, or sometimes quantitatively, accurate results at small 
computational cost. 
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Appendix A: Equivalence of Wigner transforms of products of mapping operators 

In this Appendix we show that Eqs. f )42p and f )45p are equivalent. Denoting the expression 
in Eq. PSj) for the matrix product A^y(X)B!(y' (X) by X and inserting the definition of gw(x) 
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in Eq. ({TO]) for the two g factors we obtain 



[r 



2 / + g >< r ' - 2 |mv)(mA|r + 2 )jBm(x '' X) ' (71) 
where we have used completeness on the set of mapping states. Letting r c = (r + r')/2 and 
r r = r — r', with a similar change of variables for the z variables, and using the relation 
(r — I \r' + !") = <5(r r — z c ) we find 

1 = (2^p / dTc J dpdp ' J dZcdZr Am{rc + jiP)e^ +p ' yz c /h e^- p '^ 2h (72) 
(r c + -j - 2 c |mA')(m A |r c + y + ^ c )5 m (r c - 

We have not indicated the dependence on X in this equation. Using the definition of g XX i 
we can write 

(r c + Zj t ~ z c \m y )(m x \r c + ^ + z c ) = 1 / dp c e~ i2 ^' Zcln g\x{r c + ^r,Pc)- 



4 ~ c '"^/v"-a,. c ' 2 ' c/ (2vr/i)^ 7 ^ c ' 4 

Furthermore, 

<&A'(r c + pPc) = e^/ 4 > v -5Av(r c ,p c ). (73) 
Inserting these expressions into X, integrating by parts to move the translation operator to 
the other functions in the integrand and returning to the z and z' functions, we find 

X = — J dr c dp c g xx ,(r c ,p c ) J dpdp' dzdz' ^ -^'^ 

xA m (r c + ^,p)B m (r c ~^,p'). (74) 
Next we make use of the Fourier transforms of A m and B m , 

A m (r c + Z -,v) = J dadr e^^'^+^a^a, r) 

B m (r c + ^p') = J da'dr' e^'< r ^ + ^/ n /3 m (a' , r') (75) 

Inserting these expressions into the previous form of X, performing the integrals over z and 
z' to obtain delta functions and finally performing the integrals over p and p', we obtain 

X = J dr c dp c g xx ,(r c ,p c ) [^rp J dodrda'dr' e^ + ^l h a m {a,r) 

e i ( r - CT '- T '' CT )/ 2 ftg m ( ( / ) T '^ e i((r'-rc+r'-p c )/h 



As shown in Ref. 



the quantity in square brackets is (A m B m )w, which establishes the 



equality between the expressions. 
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Appendix B: Evolution operators and projections onto the physical space 



In this Appendix we establish the equality given in Eq. ( )52|) that shows iC m commutes 
with the projection operator V. Inserting the definitions of B m (X), B^(X), p m (X) and 
pZ.{X) given in Eqs. flU, flU, & and the equation takes the form, 



dX B$ (X) [ J dx c^(x)i£ m g„, v (x)\p™ (X,t) 
dX B^ (X) / dx g^^(x)iCmC vu \x) p v ^[X,t) 



(76) 



In Ref. 



141 ] we showed that 



(77) 



J dx g^^{x)%C m c vv ,{x) = iCf, 

so that the right side of Eq. (176"j) takes the form 

dX B^'(X)tC,,^ u ,p^'(X,t). (78) 

After an integration by parts with respect to the mapping phase space coordinates, the left 
side of Eq. (176!) can be written as 



dX B$(X)i£*p^'(X,t). 



(79) 



Since = iC^^yy' (see Eq. (0J), this establishes the identity. 

Following a similar strategy we can show that the Poisson bracket mapping operator 
i£m B does not commute with V. To do this we show that 

J dX B m {X)iC p m B p v m {X , t) ± J dX BZ{X)iC™p m {X,t). 

Since iC m = iC^ B + i£' m , it suffices to show that 

J dX B m (X)i£ m pl(X,t)^ J dX Bl(X)iC' mPm (X,t), 

and we are again led to consider integrals like those in Eq. (!76j) except that iC m is replaced 
by iC' m . In Ref. (hj, Eq. (29), we established 



J dx g^{x)i£ m c vv ,{x)p v ^ '(X,t) = ^<W Tr (|| • T^r)- 
Evaluation of the corresponding integral using integration by parts gives 



(80) 



1 (dh^ L dp uv 
dx c^{x)iC! m g vlv {x)p v ^ (X,t) = -<W (-^ 7^ 



(81) 
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Thus, the evaluation of Eq. (|80|) yields 

establishing the fact that does not commute with the projection operator V . 



Appendix C: Integration scheme 

We present an integration scheme to solve the system of equations ( )58|) . This scheme is 
based on an operator-splitting method, which is motivated by the separation of time-scales 
between the electronic and nuclear motions in the problem. This method is time-reversible, 
symplectic, and includes an analytic solution for the quantum subsystem degrees of freedom. 

The formal solution to the PBME fl59|) for a dynamical variable B m (X,t) is, 

B m (X,t) = e lC ™ Bt B m (X,0), (83) 

and writing a short time decomposition of the propagator we have, e* £ ™ s * = Ylk=i e tC ™ BAt . 
The total Hamiltonian in Eq. f )36|) can be written as a sum of two parts, H m = Hi + H2, 

#i = 7^, H 2 = V (R) + ±-h xx '(R)(r x r x ,+p x p x ,). (84) 

The first part in the decomposition of H m is chosen to be the kinetic energy of the envi- 
ronment, and the second part contains the remainder of the terms in the mapping Hamilto- 
nian (|36|) . This choice of decomposition is motivated by the desire to enhance the stability of 
the approximate integration scheme and to minimize the difference between the true Hamil- 
tonian H, which is conserved by the exact dynamics, and the pseudo-Hamiltonian -ff pseu do) 
which is exactly conserved by the approximate dynamics dictated by the integration scheme. 

Partitioning the Hamiltonian in this way generates new Liouville operators, iCj = 
— \Hj, -} x , such that iCq = i(C\ + C 2 ). We then express each of the short-time propa- 
gators using the symmetric Trotter decomposition, 

el (C 1+ C 2 )At = ^C^^At^if) + 0{At 3y (g5) 

The decomposition is most useful if the action of the individual propagators e tCiAt on 
the phase points of the system can be evaluated exactly. When this is the case, the exact 
dynamics of the integration scheme is governed by the pseudo-Hamiltonian 

At 2 ( P d 2 H 2 P 1 dH 2 dH 2 \ 4 



#pseudo - H + ^ f • dRdR ■ M ^ dR ■ QR ) + 0(Af 
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For the decomposition in Eq. ( 184)) . the difference between the true Hamiltonian H and the 
pseudo-Hamiltonian is of the form of the standard Verlet scheme, and depends only on the 
smoothness of the effective potential H 2 (R) for the environment variables R and P and not 
on the smoothness of the phase space variables representing the quantum subsystem. This 
form of the splitting is particularly helpful when the system passes through regions of phase 
space where there are rapid changes in the populations of the diabatic quantum states. For 
trajectories passing through such regions, which are common in mixed quantum/classical 
systems, other decompositions of the Hamiltonian result in unstable integrators unless very 
small timesteps At are chosen. 

The evolution under iC\ gives rise to a system propagator on the environmental coordi- 
nates alone, 



JCiAt 



( r{t) \ 




( r(t) \ 


pit) 




pit) 


R(t) 




Rit) + ^At 


\ m ) 




\ P(t) J 



(86) 



Evolution under iC 2 looks somewhat more complicated; however, it may evaluated analyti- 
cally as R{t) is stationary under this portion of the dynamics. The equations of motion are 
as follows: 

h x > x '(R) dvx h x ' x '{R) 



dr x 
dt 
dP 
~dt 



(R) / dpx 
h ' dt 



h 



-r\> 



dR 
~dt 







dV (R) 1 dh x > x \R) 



(r A r A / +p x px')- 



OR 2h OR 
Consider the spectral decomposition of the mapping Hamiltonian, 



h xx '(R) = C x ,(R)E,(R)C; x \(R), 



(87) 



(88) 



where E^R) are the eigenvalues (adiabatic energies) of h. The columns of the matrix C 
correspond to the eigenvectors of h. To simplify the evolution equations for the mapping 
variables, we use the spectral decomposition of h to perform the following transformation, 



r A = C xx \r x/ 



Px 



(89) 



The two coupled equations for r A and p x from Eq. (1871) . in the tilde variables, become 



dt 



Ex{R ) : 

h 



-Px, 



dpx 
dt 



-E X (R) : 

h 



-r\. 



(90) 
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The above system may be expressed as the matrix equation, = Aiu, where the 
transpose u T of the vector u for an arbitrary quantum subsystem is written as, u T = 
(fi,Pi, • • • ,vn,Pn)- The matrix M has the simple block diagonal form, 




(91) 



where U\(R) = E\(R)/h, ® is the matrix direct sum, and a y belongs to the set of 2 x 2 
Pauli matrices. 

The general solution to Eq. fl90l) is u(t+ At) = e MAt u(t), where, in this particular case, 
the matrix exponential has the form 



^MAt 



^j^(cos(u;AAt)l + i sm(u\At)a y 



(92) 



The time evolved tilde variables are thus obtained, 

f x (t + At) = cos(u x At)r x (t) + sm(uJxAt)p x (t), 
px{t + At) = cos{uxAt)px{t) -sm{u x At)f x {t). 



(93) 



These results can then be back-transformed to the original (untilded) variables, and used to 
solve for the time-evolved momenta from equation (!H7|) . The explicit form for P(t + At) is 



P(t + At) = P{t)- d -^At 



Hence, the evolution under £ 2 is given by, 

/ r{t) \ 

e i£ 2 At PQ) 

R(t) 

\ nt) J 



( r{t + At) \ 
pit + At) 
R(t) 
V Pit + At) J 



(94) 
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Figure 1: Plots of the momentum distribution p(Pfi na i) after passage through the avoided crossing: 
QCLE (solid lines), PBME (dashed lines). The parameter values are A = 0.01, B = 1.6, C = 0.005 
and D = 1 (both panels), and the initial momentum is Po = 11 (left panel) and Po = 20 (right 
panel). All parameters are reported in atomic units. 
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Figure 2: Ground adiabatic state populations Ps t (t = 50 Js) versus 7. The quantum results are 



taken from Rcf. 



52(| and the QCLE results are from Ref. 



55J. The parameters in the FLV model 



are: uj x = 0.001, u Y = 0.00387, M x = 20000, My = 6667, a = 3, $ = 1.5, X x = 4., X 2 = X 3 = 3. 
and A = 0.01, all in atomic units. 
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— QCLE, no bath 




Figure 3: Px momentum distributions after passage through the conical intersection. The plot 
shows distributions obtained from simulations of the QCL and PBM equations for the FLV model 
without and with coupling to a bath of harmonic oscillators. The number of oscillators is Nb = 100 
and the temperature is T = 300K. 
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Figure 4: Comparison between the quantum-mechanical and full QCLE reaction reaction proba- 
bilities, with that given by the approximate PBME dynamics, as a function of the excess energy. 
Parameter values: M x = 6289, My = 8385, A = 0.00136, a = 0.458038, X = 5.0494, D e = 
0.038647, D rep = 0.02. (All quantities in atomic units.) 
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